#I think within R we should just assume we have the data in a data frame in the
# format that results from the below example. Similar to what we do in matlab
# with a dates vector and a 'data' matrix

#need included Sparkling.wtr example
d = read.table("Sparkling.wtr",sep='\t',colClasses=c(DateTime="POSIXct"),
    na.strings=c("NA","NaN"),header=T)

##############################################################################
#####  Function to downsample data to a given timestep #######################
#####  dataIn:   R data frame with time in column 1    #######################
#####  timestep: new timestep in seconds               #######################
##############################################################################
gDownsample <- function(dataIn, timestep){

	dataOut = gGapFill(dataIn,timestep)
	
  #dataOut

}

##############################################################################
#####  Function to filter data given cutoff frequencies ######################
#####  dataIn:   R data frame with time in column 1    #######################
#####  minFreq:  lower cutoff frequency (value of zero signifies no lower cutoff
#####  maxFreq:  upper cutoff frequency (value of zero signifies no upper cutoff
#####  toPlot:   boolean value (1 signifies a plot is requested)  ############
##############################################################################
gFilter <- function(dataIn, minFreq, maxFreq, toPlot){
	Ts = ts(dataIn)
	N = nrow(Ts)
	timeStep = mean(diff(Ts[1:N,1]))
	samplesPerYear = (60*60*24*365)/timeStep
	xfft = fft(Ts[1:N,2],"false")
	tt1 = seq(0,ceiling(N/2)-1,1)/N           # vector of digital frequencies
	tt2 = seq(floor(N/2),1,-1)/N
	tt = c(tt1,tt2)
	freqs = tt*samplesPerYear                 # vector of frequencies in samples/year
	lowerBound = freqs<minFreq             
	if(minFreq==0)                                         	
	{
		lowerBound = c()                               
   	}
	upperBound = freqs>maxFreq
	if(maxFreq==0)
    	{
		upperBound = c()
	}
	xfft[lowerBound] = 0                      # remove frequencies outside of given bounds
	xfft[upperBound] = 0
	data = abs(fft(xfft,"true")/length(xfft)) # create filtered vector of data
	if(toPlot)
	{
		plot(1:length(data),data,"l",c(min(1:length(data)),max(1:length(data))),c(min(data),max(data)),"y","Filtered Data","","Time","Dependent Variable")
	}
	data
}

##############################################################################
#####  Function to fill in missing data in a timeseries ######################
#####  dataIn:   R data frame with time in column 1    #######################
#####  timestep: new timestep in seconds  ####################################
##############################################################################
gGapFill = function(dataIn, timestep){

	Ts = ts(dataIn)
	newTs = seq(Ts[1,1],Ts[nrow(Ts),1],timestep)         
	span = Ts[nrow(Ts)]-Ts[1]
	numEntries = round(span/timestep)
	splined = approx(Ts[1:nrow(Ts),1],Ts[1:nrow(Ts),2],xout = newTs,"linear")
	newDates = as.POSIXlt(splined[[1]],"",origin = "1970-01-01")
	dataOut = data.frame(newDates, splined[[2]])
	dataOut
}

##############################################################################
#####  Function to remove NaNs from a timeseries #############################
#####  dataIn:   R data frame with time in column 1    #######################
##############################################################################
gRemoveNans = function(dataIn){
	NAloc = which(is.na(dataIn[2])==TRUE)
	dataIn[NAloc,2] = 0
	if(NAloc>1)
	{
		dataIn[NAloc,2] = dataIn[NAloc-1,2]
	}
	dataOut = dataIn
	dataOut
}

##############################################################################
#####  Function to compute the power spectrum of a timeseries ################
#####  dataIn:   R data frame with time in column 1    #######################
#####  toPlot:   boolean value (1 signifies a plot is requested)  ############
##############################################################################
gSpectrum = function(dataIn, toPlot){
	Ts = ts(dataIn)
	timestep = mean(diff(Ts[1:nrow(Ts),1]))
	samplesPerYear = (60*60*24*365)/timestep
	powers = abs(fft(dataIn[1:nrow(dataIn),2],"false"))^2
	N = nrow(dataIn)
	tt = seq(-N/2+.5,N/2-.5,1)/N
	scales = tt * samplesPerYear         # frequencies in cycles per year
	if(toPlot)
	{
		 plot(scales,c(powers[ceiling(N/2):length(powers)],powers[1:ceiling(N/2)-1]),"l",c(min(scales),max(scales)),c(min(powers),max(powers)),"y","Power Spectrum","","Frequency (cycles/year)","Power")
	}
	dataOut = data.frame(powers, scales)
	dataOut
}

##############################################################################
#####  Function completed in Matlab but not in R as of now ###################
#####  dataIn:   R data frame with time in column 1    #######################
##############################################################################
gTrend = function(dataIn){

}
